In [1]:
from __future__ import division
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
from itertools import product
from datetime import *
from dateutil.relativedelta import *
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
Для выполнения этого задания нам понадобятся данные о среднемесячных уровнях заработной платы в России:
In [2]:
#Reading data
wage = pd.read_csv('WAG_C_M.csv', sep=';', index_col='month', parse_dates=True, dayfirst=True)
wage.info()
In [3]:
wage.head()
Out[3]:
In [4]:
_ = plt.figure(figsize=(15,7))
_ = wage.WAG_C_M.plot()
_ = plt.title('Average nominal wage')
Проверка стационарности и STL-декомпозиция ряда:
In [5]:
_ = sm.tsa.seasonal_decompose(wage.WAG_C_M).plot()
print('Augmented Dickey-Fuller unit root test p=%f' % sm.tsa.stattools.adfuller(wage.WAG_C_M)[1])
Сделаем преобразование Бокса-Кокса для стабилизации дисперсии:
In [6]:
wage['WAG_C_M_box'], lmbda = stats.boxcox(wage.WAG_C_M)
_ = plt.figure(figsize=(15,7))
_ = wage.WAG_C_M_box.plot()
_ = plt.title(u'Transformed average nominal wage')
print('Optimal parameter of the Box-Cox power transformation: %f' % lmbda)
print('Augmented Dickey-Fuller unit root test p=%f' % sm.tsa.stattools.adfuller(wage.WAG_C_M_box)[1])
Критерий Дики-Фуллера отвергает гипотезу нестационарности, но визуально в данных виден тренд. Попробуем сезонное дифференцирование; сделаем на продифференцированном ряде STL-декомпозицию и проверим стационарность:
In [7]:
wage['WAG_C_M_box_diff'] = wage.WAG_C_M_box - wage.WAG_C_M_box.shift(12)
_ = sm.tsa.seasonal_decompose(wage.WAG_C_M_box_diff.dropna()).plot()
print('Augmented Dickey-Fuller unit root test p=%f' % sm.tsa.stattools.adfuller(wage.WAG_C_M_box_diff.dropna())[1])
Критерий Дики-Фуллера отвергает гипотезу нестационарности, НО полностью избавиться от тренда не удалось. Попробуем добавить ещё обычное дифференцирование:
In [8]:
wage['WAG_C_M_box_diff2'] = wage.WAG_C_M_box_diff - wage.WAG_C_M_box_diff.shift(1)
_ = sm.tsa.seasonal_decompose(wage.WAG_C_M_box_diff2.dropna()).plot()
print('Augmented Dickey-Fuller unit root test p=%f' % sm.tsa.stattools.adfuller(wage.WAG_C_M_box_diff2.dropna())[1])
Гипотеза нестационарности отвергается с ещё большим уровнем значимости, и визуально ряд выглядит лучше — тренда больше нет.
Посмотрим на ACF и PACF полученного ряда:
In [9]:
plt.figure(figsize=(15,10))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(wage.WAG_C_M_box_diff2.dropna()[12:].squeeze(), lags=50, ax=ax);
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(wage.WAG_C_M_box_diff2.dropna()[12:].squeeze(), lags=50, ax=ax);
Начальные приближения: Q=0, q=1, P=1, p=1.
In [10]:
ps = range(0, 2)
d=1
qs = range(0, 2)
Ps = range(0, 2)
D=1
Qs = range(0, 1)
In [11]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
parameters_list
len(parameters_list)
Out[11]:
Out[11]:
In [12]:
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in parameters_list:
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model=sm.tsa.statespace.SARIMAX(wage.WAG_C_M_box, order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [13]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())
Лучшая модель:
In [14]:
print(best_model.summary())
Её остатки:
In [15]:
_ = plt.figure(figsize=(15,12))
_ = plt.subplot(211)
_ = best_model.resid[13:].plot()
_ = plt.ylabel(u'Residuals')
_ = ax = plt.subplot(212)
_ = sm.graphics.tsa.plot_acf(best_model.resid.values.squeeze(), lags=50, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[13:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1])
Остатки несмещены (подтверждается критерием Стьюдента), стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой). Посмотрим, насколько хорошо модель описывает данные:
In [16]:
def invboxcox(y,lmbda):
if lmbda == 0:
return(np.exp(y))
else:
return(np.exp(np.log(lmbda*y+1)/lmbda))
In [17]:
wage['model'] = invboxcox(best_model.fittedvalues, lmbda)
_ = plt.figure(figsize=(15,7))
_ = wage.WAG_C_M.plot()
_ = wage.model[13:].plot(color='r')
_ = plt.title('Average nominal wage')
In [18]:
wage2 = wage[['WAG_C_M']]
date_list = [datetime.strptime("2017-07-01", "%Y-%m-%d") + relativedelta(months=x) for x in range(0,36)]
future = pd.DataFrame(index=date_list, columns=wage2.columns)
wage2 = pd.concat([wage2, future])
wage2['forecast'] = invboxcox(best_model.predict(start=294, end=329), lmbda)
_ = plt.figure(figsize=(15,7))
_ = wage2.WAG_C_M.plot()
_ = wage2.forecast.plot(color='r')
_ = plt.title('Average nominal wage')